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The turbulent boundary layer on a flat plate, with zero pressure gradient, is simulated 
numerically at four stations between Rg = 225 and Rg = 1410. The three-dimensional 
time-dependent Navier-Stokes equations are solved using a spectral method with up to 
about 10 7 grid points. Periodic spanwise and streamwise conditions are applied, and 
a multiple-scale procedure is applied to approximate the slow streamwise growth of the 
boundary layer. The flow is studied, primarily, from a statistical point of view. The 
solutions are compared with experimental results. The scaling of the mean and turbulent 
quantities with Reynolds number is examined and compared with accepted laws, and the 
significant deviations are documented. The turbulence at the highest Reynolds number is 
studied in detail. The spectra are compared with various theoretical models. Reynolds- 
stress budget data are provided for turbulence-model testing. 


1. Introduction 

The behavior of turbulent boundary layers is far from being fully understood, and accu- 
rate predictions are difficult except in the simplest cases. Scaling laws have been proposed 
and validated by measurements. They include the “law of the wall”, the “defect law” or 
“law of the wake”, and the “log law”. Other scaling laws, such as Kolmogorov’s, can apply 
but are not specific to wall-bounded flows. The scaling laws express our understanding 
of the mechanics of turbulent flows. They are also constantly used to extrapolate results 
from one Reynolds number to the other. This is necessary because the Reynolds numbers 
encountered in practical applications are much higher than can be reached in laboratory 
experiments or in direct numerical simulations. 

The theoretical foundation of the laws is often fragile; for instance several interpretations 
have been given for the “log law”. This is troublesome when one attempts to generalize the 
laws, either to more complex flows (e.g., pressure gradients) or to other quantities (e.g., 
higher statistical moments). The different interpretations, which agree in simple cases, may 
conflict in the general case. Often the available measurements are not accurate enough 
to indicate which theory is correct, and there is no consensus. Furthermore the theories 
are usually unable to predict the value of the universal constants (e. g., the Karman and 
Kolmogorov constants). For the theory of turbulent boundary layers, see in particular 


Coles (1956) and Townsend (1956, 1976). 

Another area of controversy is the behavior of the flows at low Reynolds numbers. This 
is the range of Reynolds numbers for which the flow is turbulent but does not have a very 
wide range of scales. It is turbulent in the sense of having fluctuations with significant 
energy, which is sustained by a transfer from the mean flow. The fluctuations have a 
strong efTect on the mean flow through the Reynolds stresses they generate. On the other 
hand the range of scales is not wide enough for many widely-used arguments, which are 
based on the separation of large and small scales, to apply. Typical examples are the 
inertial range in the spectrum of small-scale turbulence and the log layer near a wall; both 
properties vanish when the Reynolds number is too low. Unfortunately, all direct numerical 
simulations to date fall into the range of “low-Reynolds-number turbulence”. This makes 
the interpretation of the results delicate; it is often hard to tell whether a result is close 
to the asymptotic high-Reynolds-number value, or even whether a finite asymptotic value 
exists. One needs to distinguish between “Reynolds-number effects” and “low-Reynolds- 
number effects” . In an effort to achieve this distinction in the present study, particular care 
was taken to prevent spurious numerical effects. The simulations also cover a rather wide 
range of Reynolds numbers (a factor of 4), and the results suggest that in the Re = 1410 
flow the low- Reynolds- number effects are quite weak. For low- Reynolds- number effects 
see in particular Coles (1962), Head and Bandyopadhyay (1981), Purtell, Klebanoff and 
Buckley (1981), Murlis, Tsai and Bradshaw (1982) and Erm, Smits and Joubert (1985). 

In experimental studies, Reynolds numbers large enough to ensure scale separation are 
easily reached; for instance “long” inertial ranges and log layers have been measured. 
On the other hand some quantities, such as the pressure fluctuations within the flow, 
are very difficult to measure. Measurements on small scales are also difficult, especially 
very close to a solid wall. Thus, experimental and numerical studies complement each 
other and will both be needed in the foreseeable future. Direct and large-eddy simulations 
of turbulence were reviewed by Rogallo and Moin (1984). In the field of wall-bounded 
flows, Deardorff (1970) and Schumann (1975) studied turbulent channel flow by large- 
eddy simulation, without resolving the viscous wall region. Moin and Kim (1982) did a 
large-eddy simulation and could resolve the wall region; they obtained close agreement with 
the well-known law of the wall. Moser and Moin (1984) performed a direct simulation of 
(curved) turbulent channel flow. In contrast to large-eddy simulations, direct simulations 
do not include modeling of the eddies smaller than the grid spacing. This limits them to 
lower Reynolds numbers, but the results are thought to be more reliable, especially close to 
the wall (Moser and Moin 1984). Spalart and Leonard (1985) presented direct simulations 
of boundary layers with a range of pressure gradients. These studies did not emphasize the 
Reynolds-number dependence of the turbulent quantities and relied heavily on wall scaling 
(based on the friction velocity u T and the kinematic viscosity u) to present the results and 
compare them with experiments. 

In addition to theoretical results like the scaling laws, fundamental turbulence research 
is expected to provide quantitative information for the calibration of turbulence models. 
These models are needed for the prediction of practical flows; their present accuracy leaves 
much room for improvement. From an accurate direct simulation, one can extract all 
the quantities that are involved in a turbulence model of any complexity. This provides 
a complete test of the model and leads to suggestions about how to improve it. Here 
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one is making the assumption that the extrapolation is possible not only to different 
Reynolds numbers, but also to more complex geometries. The present study focusses 
on conventional, statistical measures of turbulence. The modern concepts of coherent 
structures in turbulence are not ignored, but the open questions regarding the behavior 
of quantities as simple as the mean velocity and the Reynolds stresses are more urgent 
both from a theoretical point of view (scaling laws) and from a practical point of view 
(turbulence models). They are also more likely to have durable answers. 

The method used for this study was applied to equilibrium boundary layers with pres- 
sure gradients by Spalart and Leonard (1985) and to sink-flow boundary layers by Spalart 
(1986b). In both studies the flows were assumed , locally, to satisfy similarity properties 
both in the wall region (law of the wall) and in the outer region (defect law). Self-similarity 
was assumed both for the mean velocity and for the Reynolds stresses. These assump- 
tions are reasonable, but are strongly justified only in the case of the sink flow. In fact 
when Spalart and Leonard compared simulations (with zero pressure gradient) at different 
Reynolds numbers, it was found that some of the scaling laws were not satisfied. Thus 
the study as a whole was inconsistent (since the global behavior of the flow negated the 
assumptions that were made locally). The deviations from the scaling laws were especially 
striking near the wall, where the terms which depend on the similarity assumptions are 
very weak. Therefore it was expected that even if these terms were altered to remove the 
inconsistency, the results would not change significantly and the deviations would remain. 

This finding motivated a generalization of the approach and a systematic, controlled 
study of the Reynolds-number effects. The new approach makes no assumption about the 
behavior of the dependent variables. The outer and inner length scales 6 and i //u r are 
still used to define a transformation in the direction normal to the wall, but as this is a 
much weaker assumption it can only have an indirect effect on the results. The procedure 
will be described in §2.3. Although it has much in common with the procedures used by 
Spalart and Leonard (1985) and Spalart (1986b) a self-contained account of the method 
will be given. 

The numerical method was described in detail by Spalart (1986a). It is fully spectral 
in space, based on Fourier series in the directions parallel to the plate and an exponential 
mapping with Jacobi polynomials in the normal, semi-infinite direction. The time integra- 
tion is second-order accurate and hybrid; it uses a Runge-Kutta scheme for the transport 
term and the Crank-Nicolson scheme for the Stokes terms. If Reynolds-number effects are 
to be studied by numerical simulation, it is essential to ensure that the different cases are 
not run with (effectively) different resolution, which could induce spurious variations. A 
similar problem can occur in experiments, for instance if a probe of fixed size is used while 
an increase in the Reynolds number decreases the scales of the turbulence. 

There are two aspects to the question of resolution. One is the size of the periods in 
the directions parallel to the wall (or equivalently the smallest wave number). It was 
decided to keep the ratio of these periods to the displacement thickness S’" the same in 
all the simulations. Thus if there is an effect of the confinement of the flow inside a finite 
period, the effect will be as independent of Reynolds number as possible. The displacement 
thickness is an appropriate macroscale of the flow and is used because, with the present 
method, it is easier to control than the boundary-layer thickness 6 or the momentum 
thickness 0. The length scale y 0 of the exponential mapping (Spalart 1986a) is also kept 
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at a constant multiple of 6\ 

The other aspect is of course the grid spacing (or equivalently the largest wave number). 
In this case the wall region is the most sensitive and it was decided that the grid spacing, 
in wall units, should be independent of Reynolds number. Therefore the number of points 
in each periodic direction is roughly proportional to the Reynolds number. In the normal 
direction the number of points was chosen in each case such that the 10** 1 point is less 
than 9 wall units away from the wall. Thus the effects of numerical truncation will be as 
independent of Reynolds number as possible. 


2. Governing equations 


2.1. Multiple-scale approximation 

The overall goal is to obtain a set of equations which, when solved with periodic condi- 
tions in the streamwise (x) direction, can provide a good approximation to the local state 
of a boundary layer that has a slow spatial development. As discussed by Spalart (1986b) 
the incentives to use periodic conditions are both numerical (the high accuracy of Fourier 
series) and physical (no need to provide turbulent inflow conditions). The idea is to use the 
fact that both the thickness of the boundary layer and the energy level of the turbulence 
vary slowly as functions of x, as illustrated in figure 1. The final product is a set of small 
“growth terms” that are added to the usual Navier-Stokes equations. 

Figure la is a sketch of the flow, with the normal direction (y) enlarged about 10 times. 
The velocity profiles at two stations illustrate the thickening of the boundary layer. The 
solid lines represent streamlines of the mean velocity field, which are at a shallow angle to 
the wall. The dashed lines represent a coordinate system that is “fitted” to the boundary 
layer. The periodic conditions will be applied along these lines. A new coordinate r\ (which 
is constant along a dashed line) replaces y as the normal coordinate; it is chosen so that 
the boundary-layer thickness and the viscous-sublayer thickness are independent of x. In 
the sink flow the coordinate lines were obvious: they were the rays converging into the 
sink. In the general case the choice is more arbitrary and will be discussed in §2.2 and 
2.3. In terms of geometry, the key quantity is the slope S of the coordinate lines. In 
the constant-pressure boundary layer, as indicated in figure la, S is positive and small, 
although larger than the slope of the streamlines. 

The dependent variables also require a change of variables before periodic conditions 
are applied. Figure lb is a sketch of a turbulent velocity component, for instance u, vs. 
x. The signal displays fast, short-scale, apparently random fluctuations as well as a slow 
variation of the mean and of the intensity of the fluctuations with x. This variation makes 
periodic conditions inadequate. However if ones writes u as the combination 

u(x,y, 2 ,f) = U (x, y) + A(x,y)u p (x,r?,z,t) (l) 

where U is the mean (over the spanwise direction 2 and time t) and the “amplitude 
function” A is proportional to the rms of the fluctuations, then the normalized signal u p 
has zero mean and its rms is independent of x. This makes periodic conditions for u p appear 
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tolerable although not fully justified. They are not fully justified because the length- and 
time scales of the signal may vary with x, so that even it,, is not a truly homogeneous 
signal. However in a boundary layer, the variation of these scales is presumed to be slow 
and to have a weaker effect than the variation of the mean and rms. 

Differentiating (1) with respect to x yields 


du 

dx 


dU A du,, 
— + A— 1 - 
dx dx 


+ 


dA 
dx v 


( 2 ) 


Since the variation of U and A with x is assumed to be slow, the first and third terms 
appear as small corrections, “slow derivatives”, denoted by U\ and u\- The second term, 
the usual or “fast” derivative du p /dx, is denoted by u x . The formulas for Ux and ux will 
be discussed in §2.3. 


2.2. Short-scale analysis 


The short-scale analysis includes the transformation of the Navier-Stokes equations from 
the Cartesian system of coordinates ( x,y,z ) to the non-Cartesian system (x,y,z), the 
inclusion of the slow derivatives, and the Taylor expansion of the resulting equations. The 
metric coefficients S and T of the coordinate transformation are defined by 






(3) 


S has a clear physical meaning: it is the slope made by the new coordinate lines with the 
wall. The meaning of rj and T is not as clear, since they depend on a normalization. In the 
end the normalization will be chosen so that, at the value of x being considered, rj and y 
coincide, giving T = 1. Until then the identity T x = S ^ will be used to express the results 
in terms of S as much as possible. 

Let (u“ ,v* ,w“) be the Cartesian velocity components. The contravariant velocity com- 
ponents (u,v,u>) associated with ( x,y,z ) are used to preserve the form of the transport 
terms, that is to have u*d x + v”d y + w~d z = ud x + vd r , +wd z . They are defined by w = w* 

and , „ ... 




The continuity condition (it* + v’ y + =0) becomes 


(4) 


S T 

+ yi + y i = o. 


(5a) 


The x-momentum equation becomes 


ut + UU X + VU V + WU Z = -p x + — p v + 


v I W xx + 


i + s 2 „ 

~rp2 


, 2SS " + 

j, U xr , + f T ' T 2 y3 


) Hf) "h 


(5b) 


where p is the kinematic pressure. Similar terms enter the other equations. Now that all 
the derivatives have been taken, the normalization can be applied, which will simplify the 
equations: rj is identified to y and T is set to 1. The equations become 
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ut + uux + tm y + ttm z = 


u x + Vy + w z + SyU = 0, (6a) 

-p x + Sp y + l'(u xX + (l + S‘ 2 )Uyy-2SU xy + (SSy- S X + SSy)U y +U zz ). 

(6b) 

Let us proceed with the Taylor expansion of the equations. The velocity and pressure are 
split into mean and fluctuating quantities: e. g u = U + u, and the fast and slow deriva- 
tives defined by (2) are introduced. The friction velocity u T /Uoo is a small quantity and 
tends to 0 (although slowly) as the Reynolds number tends to oo. Let u T /U <*, be of order t . 
The velocity fluctuations u,v,w are of the same order. The pressure fluctuations p/U ^ are 
of order c 2 . The momentum equation, dO/dx = u 2 /£/£>, shows that the rate of change of the 
boundary-layer thickness 6 in the x direction is of order £ 2 . Then S and 6 Ax /A are of order 
e 2 . All of the slow derivatives such as U\ and u\ are two orders of magnitude smaller than 
the quantity itself, whereas the fast derivatives are of the same order as the quantity (tak- 
ing 6 as the length scale). For instance u/U Q 0 = 0(e), u x 6/U ^ — O(f), = 0(e 3 ). 

The Taylor expansion of (6) up to order e 3 is as follows (note that in each equation 
the leading terms are of order £, so that the relative order of the corrections is £ 2 ). 


V y + [U x + S y U] = 0, (7a) 

u x + Vy + w z + [u x + s y u] = o, (7b) 

u t + (U + u)u x + v(U + u) y + wu z + [UUx + VU y + Uux + U x u + ku y ] = 

- P* ~ [Px] + HUyy + V 2 «), (7 C ) 

v t + (U + u)v x + Wy + wv z + [Uvx + (Vy + 2 U S y )v + V v y ] = -p y + uV 2 v, (7d) 

wt + {U + u)w x + vw y + ww z + [Uwx + Vwy] = -p z + vV 2 w. (7e) 


The corrections in the viscous terms have been neglected. The mean and fluctuating 
components of the continuity equation are separated for clarity (7a, 7b) and all three com- 
ponents of the momentum equation are shown. These equations are the Navier-Stokes 
equations, as usually written for a shear flow that is homogeneous in x and z, with the 
addition of “growth terms” that are indicated by square brackets and can be interpreted 
as follows. The term Ux + S y U in the continuity equation (7a) is a consequence of the 
streamwise evolution of the mean-velocity profile. In the sink flow, it was identically zero 
(Spalart 1986b). In the flat-plate flow, it is zero in the wall layer and positive farther from 
the wall, so that V becomes negative: nonturbulent fluid is entrained into the boundary 
layer (in the new coordinates the dashed lines of figure la are parallel and the solid lines 
are directed towards the wall). The meaning of the ux + S y u term in (7b) is not as 
clear (it was zero in the sink flow as well). Actually the numerical method used, which 
employs basis functions that satisfy u x + v y + w z = 0 (causing the pressure term to be 
eliminated), does not allow the ux + S y u term to be included. Thus the expansion is not 
fully third-order accurate in £, but this effect is not thought to be very significant (the 
integrated contribution of the pressure term to the turbulence energy is of order £ 5 while 
the production term, for instance, is of order c 2 ). In the momentum equation, the terms 
UUx + VUy and P x are mean momentum-transport and pressure terms. The term Px is 
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independent of y, to order e 3 , and satisfies Bernoulli’s equation in the freestream: for large 
y, UU X = ~Px (in the present application, Px = 0). The terms Uux +Vu y , Uvx + Vv y , 
and Uwx + V w y represent transport, or advection, of the turbulence by the mean flow. 
The term 2U S y v in (7d), which is due to the expansion of the coordinate lines, seems 
harder to interpret. 

Finally, the terms U xu and V y v arise from the straining of the turbulence by the mean 
flow. However it would be misleading to imply that straining is accurately represented, 
because the numerical domain does not become deformed in time as in, for instance, Ro- 
gallo’s homogeneous-turbulencesimulations (1981). This weakness of the present approach 
is, ultimately, due to the fact that the approximation in (1) was based purely on energy 
considerations, without consideration of structural aspects. If one examines the influence 
of S , one finds that the mean growth term UUx + VU y is not sensitive to S , because of 
the definition of the contravariant velocities u and v (if S changes, the quantities UUx 
and VU y change, but their sum remains equal to U’U~ +V~U~). For the fluctuations 
the situation is more subtle. It can be shown that the global effect on a Reynolds stress 
such as < u 2 > behaves like the effect on the mean flow, that is, a change in S causes 
only a transfer between the terms U < u 2 >x and V < u 2 > y . However, locally the value 
of S does matter because «x and u y do not have the same phase. It is not obvious that 
the equations (7(a-e)), which were derived entirely using Reynolds-averaged concepts, will 
have the right effect on the flow structures. The production of turbulent energy due to the 
mean shear U y is well represented (it does not depend on the spatial development). The 
entrainment process by which irrotational fluid is absorbed into the turbulent region seems 
to be well represented, even locally. On the other hand, it is unlikely that flows in which 
a straining associated with the spatial development is dominant would be simulated very 
accurately. Flows in which length or time scales evolve rapidly in the x direction should 
also be avoided, as mentioned earlier. 

2.3. Long-scale analysis 

The long-scale analysis provides the values of the quantities S, Ux, «x, vx, and wx that 
enter the growth terms in (7). In the method used for sink-flow boundary layers (Spalart 
1986b) the long-scale analysis was simple; the flow was assumed to be self-similar both 
in terms of its length scales (all proportional to the distance Xq from the sink) and its 
velocity scales (all proportional to the edge velocity Q/X o). Thus one had the equations 
S — -y/X o, Ux = U/X o, ux = u/X 0 , and so on. 

In more general flows one has to consider at least two length scales: the wall length scale 
v /u T and the boundary-layer thickness 6. There are also two velocity scales: the friction 
velocity u T and the edge velocity U^. Spalart and Leonard’s (1985) study of boundary 
layers with pressure gradients was restricted to “equilibrium” boundary layers in the sense 
used by Clauser (1954), in which the outer part of the velocity profile, expressed as a 
defect and normalized by u r , is a function only of y/6. This is the “defect law”. The wall 
region was also assumed to display the universal property known as the “law of the wall”. 
Thus the natural coordinate r\ was (a function of) y + in the wall region and (a function 
of) yjS in the outer region. A smooth blending between the two regions was made. The 
fluctuations were assumed to scale with u T both in the wall region (for constant y + ) and in 
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the outer region (for constant y/6). This assumption was consistent with the assumption 
made for the mean velocity and with the literature. 

Spalart and Leonard’s analysis has several shortcomings. It is limited to “equilibrium” 
boundary layers; such boundary layers are not only rare in practice, they are also not 
rigorously defined. The “constant /?” definition (Clauser 1954) is not fully justified. It 
was also shown by Coles (1962) that even in the simplest and best explored of equilibrium 
boundary layers, the constant-pressure flow, the defect law is not satisfied at low Reynolds 
numbers, such as the one used for the simulations. Also, the overlap of a law of the wall 
and an outer-layer law requires a log layer in the mean velocity and a constant layer in the 
Reynolds stresses. The log layer was obtained by Spalart and Leonard, but the constant- 
stress layer was not. Although this could be a low-Reynolds-number effect, it is hard 
to explain why the log layer was well indicated at the same Reynolds number. Finally, 
simulations conducted at different Reynolds numbers yielded results in which the Reynolds 
stresses, near the wall, did not scale with u T and y + . Instead they showed a consistent 
tendency to rise with the Reynolds number. 

These various shortcomings indicated that (except for the sink flow) a valid study of 
boundary layers at low Reynolds numbers or with general pressure gradients should do 
away, at least, with the assumptions of a defect law and of the wall scaling of the Reynolds 
stresses. This led to the idea of obtaining the information needed to prescribe 5, Ux, 
and so on, not from some assumptions but directly by conducting simulations at several 
stations of the same boundary layer. The X-derivatives are obtained by taking differences 
between the different stations. This approach is more general, and is conceptually simpler. 
On the other hand it is more expensive, since several simulations have to be conducted. 
Also, the simulation of any one station requires upstream information. Thus one loses one 
of the advantages of the original procedure; however, the upstream information needed is 
reduced to mean-velocity and Reynolds-stress profiles. Time-dependent turbulent inflow 
values are not needed. 

The first task is to define the coordinate lines. The quantities y + , near the wall, and 
y/6, away from the wall, appear as natural coordinates. The idea is to minimize the 
nonhomogeneity along a constant-77 line, so that the growth terms are as small as possible. 
The new coordinate 77 is defined explicitly as a weighted average: 

_ y2( 10_3 y + ) + y p {y/ 6 ) 
vl + y 

where yf = 15, y 3 /6 = 0.3, y 2 = -y/ylyi and p = 5/ log 10 (y 3 /j/i). This rather arbitrary 
definition was chosen to ensure that p is a monotonic function of y, is equal to 10~ 3 y + for 
small y and to y/6 for large y, and makes a smooth transition between the two. 

Now that r? is defined, if one knows the state of the flow at two values of X, e. g., Xj 
and X 2 , one computes an approximation to 5 = dy/dX by taking differences: 

s * (9) 

x 2 -x t 

Typically the flow at the upstream station, X l5 is known from a previous simulation and the 
flow at the downstream station, X 2 , is being computed. Similarly Ux is approximated by 


(8) 


8 



\U{X 2 ,ri) - U(X u t])}/ (X 2 - Xi). For u, A x is given by 

[ u rms (X 2 ,ri)/u rms (Xi,T 7 ) — l]/(X 2 — Xi), and the formula for w is the same. Finally, 
the formula for v is slightly different: A x « [v rms (X 2 , ri)/v rms (Xi, rj) - l]/(X 2 - Xi) - S y . 
The — S y term is introduced by the normalization of v in (4). This cancels half of the 
2US y v term of (7d). The only information one needs at Xi is U , u rrns , v rms and w rms as 
a function of y. During the transient regime, one can use the values of U , u rms , and so 
on, at X 2 at one time step to evaluate the growth terms for the next step. No stability 
problems were encountered. 


3. Results 


3.1. Choice of the parameters 

The periods in the z and 2 directions are A x = 1006’ and A 2 = 256’ . The spacing 
between collocation points, in wall units, is Ai + « 20, with A 2 + « 6.7, and may vary by 
±5% from one case to another. In the y direction, as mentioned earlier, the number of 
points is adjusted so that there are 10 (nonuniformly spaced) points within 9 wall units of 
the wall. As a result of these requirements and of the need to factorize the length of the fast 
Fourier transforms, the number of points in the z, y, and 2 directions is (128 x 50 x 96) 
at R s * = 500, (256 x 64 x 192) at R s - = 1000 and (432 x 80 x 320) at R s - = 2000 . 
This last case has about 1.1 x 10 7 grid points, so 3.3 x 10 6 modes, (using the 2/3 rule 
in each direction, Spalart 1986a). The time step is adjusted so that the maximum local 
Courant-Friedrichs-Levy (CFL) number (see Spalart 1986a for the exact definition) is 2. 
As a result the nondimensional time step AtUoo/S * is about 0.3, 0.14, and 0.07 at Rg~ = 
500, 1000, and 2000 respectively. In wall units, the nondimensional time step Xtufjv is 
about 0.43, 0.32, and 0.29 respectively. The time sample for the statistics is of the order of 
2006’ /Uoo. Larger samples would help improve the smoothness of the spectra, but would 
be very expensive to generate. 

These values were chosen by monitoring the spectra, both at their lower end (to choose 
the period) and at their upper end (to choose the grid spacing). Alternatively, the two- 
point correlations can be used (Spalart 1986a). Two tests were also performed for the 
(relatively inexpensive) case Rg - = 500 to further validate the choice of the parameters. 
In the first test, the values of A x and X z were doubled, then halved, and the grid spacing 
left unchanged (thus the number of points in x and 2 also changed by a factor of 2). In 
the second test, the simulation was rerun with both much finer, and then much coarser 
resolution (the number of points in z and 2 was doubled, then halved; the number of points 
in y was varied from its basic value of 50 to 64, then 40). Since the time step is always 
regulated by the value of the CFL number, it is also reduced when finer resolution is used, 
so that both the spatial and temporal errors are reduced. 

Figure 2 displays the sensitivity of the mean-velocity profile to the numerical parameters. 
The difference between the basic simulation and the “improved” ones (larger periods or 
finer resolution) is small. This figure allows one to estimate the remaining numerical 
uncertainty. On the other hand, the simulation with coarse resolution gives significantly 
different results, including a much higher value of Ut/Uoo . Thus the resolution chosen 
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appears to be adequate, but not wasteful. Finally, the simulation with reduced periods 
gives results very close to the basic simulation, which indicats that the periods that were 
chosen to obtain satisfactory long-range behavior (low correlations at a distance A/2) 
are more than sufficient to obtain a satisfactory velocity profile. The sensitivity of other 
statistical quantities up to third-order moments was also found to be small. For fourth- 
order moments, the u and w components showed a moderate sensitivity, but the flatness 
factor of v was very sensitive near the wall: it varied from about 4 with coarse resolution 
to over 30 with fine resolution. Away from the wall, the profiles agreed well again, even 
the v-flatness. Thus one should consider statistics for moments beyond third order with 
some caution, at least near the wall. 

3.2. Streamwise evolution of the mean flow and the shear stress 

Four stations of the boundary layer were simulated, with Re- = 400, 500, 1000, and 2000 
respectively. The values of Re are approximately 225, 300, 670, and 1410. Detailed results 
from the first station, Re- = 400, will not be shown. It was computed using the approach 
of Spalart, and Leonard (1985). While this approach was shown to have some deficiencies, 
it is satisfactory for the generation of inflow data. A “normal” turbulent boundary layer 
was obtained except for exaggerated low-Reynolds-number effects. 

The possibility of sustaining turbulence at Re — 225 is in disagreement with Preston’s 
(1957) estimate that the lowest Reynolds number at which “fully developed turbulent flow” 
can occur is 320. It is based on a comparison between pipe-flow and boundary-layer data, 
and on the rather abstract argument that when the length of the log layer is reduced 
to zero (because the inner and outer regions overlap), fully developed turbulence cannot 
exist. One may need to distinguish between the concepts of “fully developed” turbulence 
and of “sustained” turbulence. Indeed, in the present results for Re below about 600, a 
normal log layer is not observed (see below). However the turbulence was sustained. A 
comparison with sink-flow boundary layers may be useful. These flows were also studied 
with the present method (Spalart 1986b). The threshold Reynolds number was found to be 
Re = 330, which is in good agreement with experimental data. In sink-flow simulations, 
the turbulence quickly collapsed when the Reynolds number was lowered to below 330; 
thus the present method seems reliable. The sink flow has a strong favorable pressure 
gradient, which has a stabilizing effect (it is known to inhibit transition and to induce 
relaminarization). When this stabilizing effect is removed, turbulence should be sustained 
at lower Reynolds numbers. Curiously, Preston’s argument about the overlap of the inner 
and outer region predicts the opposite effect: the limit would be lower with a favorable 
pressure gradient. Thus Preston’s estimate for the sink flow would be significantly lower 
than the value of 330, which has now been obtained by two independent approaches. 

Figure 3 shows the growth of the displacement- and momentum thicknesses in the 
streamwise direction. Since the origin is arbitrary in the x direction, the notation A R x is 
used to represent the Reynolds number based on the distance from the inflow station and 
on Uoo. The scale is enlarged in the normal direction. The momentum-balance equation, 
d(Re)/d(R x ) = c//2, is checked by drawing segments centered on the Re points with a 
slope of Cfj 2. They show that the balance equation is well satisfied and that the spacing 
between stations is small enough. In figure 4 the dependence of the friction coefficient cj 
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and the shape factor H = 6~ /Q on R$ is shown; it compares well with experimental data 
(Coles 1962, Purtell et al. 1981, Murlis et al. 1982, Erm et al. 1985). However as R$ 
increases, the cj decreases slightly less than the experiments show. Figure 2 suggests that 
the uncertainty on U^/ur is of the order of ±0.25, or ±1.25%. This translates into ±2.5%, 
or ±0.12, for 10 3 c/, which is of the same order as the differences in figure 4a (for R$ = 670 
and larger). The Re — 1410 simulation was continued with slightly improved resolution 
(96 points in y instead of 80); the cj showed no tendency to decrease. The error could also 
be caused by effects neglected in the multiple-scale approximation; unfortunately it does 
not seem possible to estimate these effects quantitatively. 

The velocity profiles are plotted in figure 5a using wall variables. Experimental results at 
Re — 617 (Erm et al. 1985) and 1368 (Murlis et al. 1982) are also plotted and may be 
compared with the present results at R$ = 670 and 1410. For Re equal to 670 and larger, 
the profiles closely follow the logarithmic law [ 7 + = log(y + ) / k + C , with constants k = 0.41 
and C = 5, starting at y + rs 30. For Re — 300, the curve is significantly higher in that 
region; however it is still rather straight. This raises the question of exactly how the log 
layer is defined. 

A logarithmic layer is a region in which the quantity dU /d{\og y), or ydU /dy, is constant 
and equal to u t /k. In high-Reynolds-number boundary layers, this region is long and 
bounded on both sides by regions in which ydU/dy takes values larger than the constant. 
Thus the logarithmic layer and the value of u t /k can be found by seeking the minimum of 
ydU/dy vs y (this amounts to seeking the inflection point in figure 5a). The value of C can 
then be determined at the same position. This seems to be the only rigorous definition. 
The usual procedure of adjusting k and C until the straight line fits the velocity profile 
“well” in the coordinates of figure 5a is not exempt from arbitrariness. Of course, the 
procedure based on ydU/dy is more sensitive to noise, and will yield an erroneously low 
value of u t /k for noisy data (C may then be overpredicted). Thus it may not be applicable 
to experimental data. However, in numerical results, the mean-velocity profile shows a low 
level of noise once the sample is sufficient. The profiles of y + dU + /dy + are shown in figure 
5b and are seen to be smooth enough in the region of interest. Note that if the value of k 
is assumed to be known, one can deduce the value of u T (and hence the wall stress) using 
only measurements in the log layer. This is essentially the method of the “Clauser plot” 
(Clauser 1954). 

When the Reynolds number is low, the quantity ydU/dy still exhibits a local minimum, 
but the minimum is narrow (figure 5b). Several interpretations are possible. The first is 
an “infinitely short” log layer. Figure 5b shows that at Re = 300, this log layer has a 
low value of about 0.315 for the “apparent Karman constant” k' . This is reminiscent of 
Simpson’s interpretation (1970), although Simpson made deductions from measurements 
taken in the outer region of the flow. He proposed the law k' = K 0 (Re/6 000) ^ with k 0 
the asymptotic value, about 0.40. This yields 0.28, 0.31, and 0.34 for our three cases. 
The present results indicate values of 0.315, 0.380, and 0.404 for k' . They show a smaller 
departure of k' from k 0 , and the value of k? is essentially unaffected down to Re « 670. 
The apparent values C', in the numerical results, are 2.51, 4.25 and 4.74 for the three 
cases. The uncertainty is estimated to be about ±0.01 for /c' and ±0.25 for C' . 

The second interpretation is simply that the log layer exists at Re = 670 but has dis- 
appeared at Re = 300. The log layer normally covers the region between y + « 30 and 


11 


y / 6 « 0.15. These two positions coincide when rs 200, which corresponds to Re « 400. 
The wall and wake regions begin to overlap for Re lower than 400, so that the local mini- 
mum of y 4 dU + /dy + is reduced to a point and is higher than its normal value. The results 
shown in figure 5 imply that at Re = 300, the Clauser-plot method yields an erroneously 
high “apparent friction velocity” u' T if the decrease in k' is ignored ( u' T equals u t ko/k' 
in which u T is the true friction velocity). Curiously, the “apparent log layer” obtained 
by fitting a straight line through the inflection point of the profile (y + ss 40) follows the 
velocity profile closely up to a high value of y/ 6 : about 0.35 as opposed to 0.15 normally. 
Purtell et al. (1981) also remarked on the tendency of the log layer to extend to larger 
values of y /6 at low Reynolds numbers. 

The “strength of the wake” A U + is defined by Coles (1962) as the maximum value of 
([/ + — £7+ ), the deviation of the velocity profile over the log law (which occurs in the 
region near y /6 — 0.8). The quantity is constant, with a value of about 3, at high Reynolds 
numbers when the defect law is satisfied, but decreases for values of Re lower than about 
6000 (Coles 1962). This is a prime example of low-Reynolds-number effect. In defining 
A U + at low Reynolds number one has to choose between the “apparent” log layer and the 
normal log layer (/c = 0.41, C — 5). Furthermore the value of A U + is very sensitive to the 
constants used for the log layer. If C is changed from 5.0 (Coles 1962) to 5.2 (Murlis et al. 
1982) A U + is clearly reduced by 0.2; if k is changed from 0.40 to 0.41 A t/+ is increased by 
about 0.4 (at Rq — 1410) which is about 25%. Figure 6 compares the computed values of 
A U + , using both definitions, with Coles’ curve and other experimental results by Murlis 
et al. (1982) and Erm et al. (1985). The curve based on the apparent log law rises like 
the experimental curves, but the trend in the other curve (based on k. — 0.41, C = 5.0) is 
very different: it is essentially flat with values of about 1.4, much lower than the accepted 
value of 3. This trend was already apparent in figure 4a. Very accurate measurements or 
simulations over a wide Reynolds-number range, as well as a strong consensus on the value 
of /c, will be needed before definitive results are obtained for A U + . 

A reliable definition of the boundary-layer thickness 6 is needed to express the results 
in the upper part of the flow. The definition of 6 as the point where U/U^ takes a given 
value (typically 0.99 or 0.995) is not very satisfactory. It depends on small differences; 
it is inconsistent with the defect law since u T depends on x. The situation is especially 
confusing at low Reynolds numbers since the velocity profile is not invariant, whether it 
is normalized by or u T . On the other hand, the shear-stress profile is very close to 
invariant when normalized by its wall value as will be shown. Let r(y) be the total stress 
and r(y)' 4 " denote r(y)/r(0); t + varies smoothly from 1 to 0 as y varies from 0 to oo. Let 
the “stress thicknesses” dj and 62 be defined by 

r 00 r 00 

61 = r + (y)dy and 6 2 = r + (y) (l - r + [y))dy. (10) 

Jo Jo 

These definitions were made by analogy with the definition of the velocity thicknesses 6 * 
and 6 . One can define a “stress shape factor” H T = 61 / 62 . For reference, the shape factors 
of a triangular and a rectangular stress distribution are 3 and Too, respectively. The 
Gaussian exp(-y 2 ) and the cubic (1 — 3y 2 T 2y 3 ), although both bell-shaped, give quite 
different values: 3.41 and 3.88. Thus H T is a rather sensitive measure of the shape of the 
stress distribution. The values of H T at the four stations of the boundary layer are about 
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4.0, 3.9, 4.0 and 3.9; they are close to each other, and to the cubic. In contrast the velocity 
shape factors, H and Clauser’s shape factor G = (Uoc /u r )(H — 1 )/H (which was designed 
to be invariant), are much more sensitive to low-Reynolds-number effects. See figure 4b 
for the values of H; G varies from 7.4 to 6.5. This suggests using the stress, instead of the 
velocity, to define the thickness 6. In order to relate this concept to the familiar concept 
of 6 as the edge of the turbulent region, or of the region with mean shear, the following 
definition was adopted: 

6 = 1.85 6,. (11) 

The factor 1.85 gives a good match with published values, for instance Klebanoff’s (1954). 
This definition has the advantage of being an integral instead of a local quantity; on the 
other hand, it applies only to the zero-pressure-gradient boundary layer. Figure 7 shows 
the shear-stress profiles plotted versus y/6. The agreement with Klebanoff and the collapse 
of the total-stress profiles are very good, thus confirming the behavior of H r and justifying 
the definition of 6. At y = S the velocity U /Uoc is equal to 0.9965, 0.9974, and 0.9977 for 
Re = 300, 670 and 1410 respectively. The corresponding values of the defect (/Too — U)/u T 
are 0.065, 0.052, and 0.050. Finally, the values of r + at y = 6 are 0.028, 0.026, and 
0.023. Definitions of 6 based on any of these three “reasonable” candidates were tried and 
produced values that differ from (11) by up to 5% and significantly degrade the collapse 
of the curves in figure 7, especially between 0.56 and 0.86. 

In figure 7 the total-stress profiles have zero slope at the wall, as they should since the 
pressure gradient is zero, and one would naturally expect a parabolic behavior for t + near 
the wall. However the profiles, at least for Re = 670 and 1410, show a flat part around 
0.16 with finite slope, dr + /d(y/ 6) of about -0.5. A simple argument is outlined in the 
appendix, which agrees with Townsend’s 1956 analysis. It suggests that this behavior is 
indeed correct, and that at high Reynolds numbers the total-stress profile approaches the 
wall with a finite slope of the order of -0.6, with the slope falling to zero only within 
the buffer layer (which becomes very thin compared with 6). This finding is significant, 
because it means that even the zero-pressure-gradient boundary layer does not contain 
a constant-stress layer. The stress behaves in a manner similar to a favorable-pressure- 
gradient flow, and the value -0.6 is of the same order as the value in a channel flow, namely 
-1 (defining 6 as the channel half-width). Again, as in the sink flow, a layer of constant 
stress and one of logarithmic velocity profile are incompatible (Spalart 1986b). 

Figure 8 shows the velocity profiles, now plotted versus y/6. In figure 8a the velocity U 
itself is plotted and normalized by U^. Again, experimental results at Re =617 and 1368 
are plotted, and the agreement is quite good (for Erm et a/.’s flow the following values were 
assumed: it r /^oo = 0.0505, 6 + = 324; for Murlis et al U t/Uoo = 0.0442, 6 + = 547). In 
figure 8b the velocity defect ( U — U^/ur is plotted; in addition to Erm et a/.’s and Murlis 
et a/.’s results, the experimental results of Klebanoff (1954) at high Reynolds numbers are 
shown (assuming u T /Uoo = 0.0375). Again the agreement is acceptable but there is no 
trend towards Klebanoff’s high-Reynolds-number curve, as Coles’ theory would predict. 
The numerical results seem to satisfy a defect law, but one that does not quite agree with 
Klebanoff’s. Note that the results as presented in figure 8 depend on 6, which is not firmly 
defined. Note also that using the apparent friction velocity u' T (deduced from a Clauser 
plot) instead of the true u T to normalize {U — U^) in figure 8b would create a trend towards 
Klebanoff’s curve and thus improve the agreement with Coles’ model. 
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3.3. Streamwise evolution of the turbulence 


The behavior of the Reynolds shear stress - < uv > is predictable. As shown in figure 
7 the normalized total stress r + , plotted versus y/6, varies very little with the Reynolds 
number. In the same figure, the Reynolds shear stress — < ur > + shows a weak sensitivity 
to the Reynolds number for y + larger than about 50. This dependence is primarily due 
to the variation of the (small) viscous stress. Near the wall, of course, the Reynolds 
stress falls to 0 at a different value of y/6 depending on the sublayer thickness. Figure 
9 shows both stresses and the turbulent-energy production - < uv > + dU + /dy + in the 
wall region. The stresses are rather insensitive to Reynolds number up to a y + of about 
15, where the low-Reynolds-number profiles bend down. The outer region, in which the 
stress falls to zero, impinges on the inner region, in which the total stress is essentially 
constant. The higher-Reynolds-number Reynolds-stress profile agrees very well with Kim, 
Kline and Reynolds’ (1968) experimental curve. The collapse of the production profiles 
is striking. They agree with each other, with Kim et o/.’s curve, and with the theoretical 
profile 1/k y + even in the region where the Reynolds-stress profiles are far from collapsing. 
At low Reynolds numbers the decrease of — < uv > + and the increase of dU + /dy + (figure 
5) cancel each other to a remarkable degree. A simple explanation of this fact has not 
been found. 

The behavior of the normal Reynolds stresses < u 2 >, < v 2 > and < tv 2 > is more 
complex than that of the shear stress. The rms values of these fluctuations are plotted 
in figure 10 vs y/6 , nondimensionalized by u T , and compared with Klebanoff’s (1954) 
results. The computed curves agree with Klebanoff’s curves reasonably well except in 
the wall region, below t/ + = 50. Outside this region the curves, at least for v and w, 
show a moderate but consistent rise with Reynolds number. Perry, Lim and Henbest 
(1985, see also Perry, Henbest and Chong 1986) explain this rise by the lengthening of 
the inertial range towards higher wave numbers, the low-wave-number part of the spectra 
being independent of Reynolds number. This hypothesis is supported by figure 11. Power 
spectra of the three velocity components are shown at y — 6/2, nondimensionalized by u T 
and 6, for the Re = 300 and 1410 cases. Spanwise spectra are shown, as they tend to be 
smoother than the streamwise spectra. Although the sample is marginal for the lowest 
wave numbers, the indication is that with a larger sample the spectra will coincide well up 
to a k z 6 of about 10. 

The theory implies that by extending the inertial range to infinite wave numbers one 
should obtain the “infinite Reynolds number” value for each Reynolds stress at each value 
of y/6. Such an extension would add to the stress a correction Cj(ufc/u r ) 2 in which C 1 is 
a universal constant, u*. is the Kolmogorov velocity scale (ut ) 1 / 4 , and t is the dissipation 
rate. A value of 2.8 was assumed for Cj, which is smaller than Perry’s value. This will 
be discussed later with the spectra (§3.4). For the computed flows the dissipation rates 
were available, whereas for Klebanoff’s results the dissipation was estimated by scaling the 
dissipation computed in the Re = 1410 simulation. Figure 12 shows that the correction 
dramatically improves the collapse of the profiles except for < u 2 > at Re = 300, and near 
the boundary-layer edge (y/6 > 0.8). In that region the turbulence is intermittent; the 
irrotational fluctuations, which do not follow Kolmogorov’s theory, contribute a significant 
part of the energy. The agreement with Klebanoff’s corrected values is also rather good. 
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The profiles still do not collapse near the wall, for y + less than about 50, which could be 
expected since no inertial range exists in that region. The obvious similarities between the 
deviations observed near the wall suggest a model of the type f\(y/6) + / 2 (y + ) where / 2 
is zero for y + larger than about 50; however there is little theoretical support for such a 
model (one could also try a product, fi{y/S) x / 2 (y + ))- 

The results in figure 12 confirm the validity of the scaling based on 6 and u T , and of 
Perry et al.'s viscous correction (although with different constants). The possibility of 
extrapolating moderate-Reynolds-number simulation results to higher Reynolds numbers 
using just a simple argument and one universal constant C 1 is very attractive (however 
the correction still does not apply in the wall region, and the theory still cannot predict 
the value of cy). 

When shown in wall variables the profiles also rise with Reynolds number (figure 13). 
The behavior of w is the most striking: the rise is significant even very near the wall. Ex- 
perimentally, Coles (1978) and Purtell et at. (1981) observed that the rms of u increased 
with Reynolds number for y + larger than 15. Coles also reports a wide scatter in the w 
values, but without a clear trend versus Reynolds number (1978, and personnal communi- 
cation 1985). The results in figure 13 are in sharp disagreement with the concept of wall 
scaling (based on u T and is), a concept that was very successful when applied to the mean 
velocity. By directly extending the arguments from the mean velocity to the stresses one 
would predict a good collapse of all curves in figure 13, especially near the wall. Note 
also that a viscous correction by Cj(ufc/u r ) 2 could not improve the collapse in figure 13 
because, at a given y + , Uk/u T varies very little with Reynolds number. 

Figure 14 shows the rms of the pressure fluctuations, normalized by u 2 . The behavior 
is very similar to the velocity rms: a slight rise with Reynolds number in the outer region 
but a strong rise near the wall. The wall value rises from about 1.9 to 2.7, which is well 
within the experimental scatter. This trend is in agreement with experimental results 
reviewed by Willmarth (1975), and is consistent with Townsend’s (1976) law for the wall 
value, < p 2 > ju A = Clog(6 + ) + D\ the comparison between the cases Re = 670 and 1410 
suggest C' « 2.05. The vorticity intensities collapse well when plotted versus y/6 and 
normalized with the “mixed” scale yju^/dis (which is appropriate for the dissipative scales 
of motion), but again show a significant rise with Reynolds number when plotted in wall 
units. This lack of collapse of the velocity, pressure, and vorticity intensities near the wall 
has important implications for theories and turbulence models. Apparently, one cannot 
rely on a “law of the wall” for these quantities as for the mean velocity. Note that the 
multiple-scale approximation that was made is probably most valid near the wall, where 
the natural length- and time scales are the smallest. The growth terms are also very small 
near the wall. Thus, it is unlikely that the surprising behavior of the turbulence statistics 
near the wall could be due to the approximations that were made. 

The Reynolds-number effects in figures 13 and 14 can be interpreted in terms of the 
theory of “active” and “inactive” motions proposed by Townsend (1961) and Bradshaw 
(1967). This theory was initially developed to explain the observation that some of the 
triple correlations were Reynolds-number dependent near the wall. The active motion 
“produces the shear stress and its statistical properties are universal functions of r and 
y” ; the inactive motion is “effectively irrotational” and does not produce shear stress 
(Bradshaw 1967). The inactive motion has length and time scales that are large compared 


15 



with the viscous-layer scales, and is not directly connected with the wall shear stress. 
Therefore it may not scale with u r and u (Townsend 1976, p. 138). 

Figure 15 shows the length scales responsible for the extra energy at higher Reynolds 
numbers. The power spectra of the four nonzero components of the Reynolds-stress tensor 
are shown. The wave number and the energy are nondimensionalized using wall variables 
(u T and i/). Spanwise spectra are presented, in the y + = 15 plane, for the cases Re = 
300 and 1410. Figure 15 shows that as the Reynolds number increases, the energy of 
fluctuations with wavelengths less than about 150 wall units is not affected. In the x 
direction, only waves longer than about 300 wall units show deviations. Thus, wall scaling 
is violated only by those scales of motion much larger than the thickness of the wall layer 
(roughly 50 wall units). These large-scale motions contribute most to < w 2 > and little 
to < v 2 > and — < uv >. This agrees well with the concept of inactive motion. 

A simple model of the inactive motion is that the pressure disturbance following a 
large coherent structure, combined with the no-slip condition, create a situation similar 
to Stokes’ oscillating boundary layer (Schlichting 1979). From the known solution to 
this equation one can predict the dependence of the extra energy on y + . If one takes 
a typical wave-number « 0.005, for which the spectra exhibit inactive motion, and 
a typical convection velocity for a large structure c + « 15 (or 0.75[/oo), one obtains a 
frequency n + zz 0.075. In figure 16, the rms profiles corresponding to the Stokes solution 
at n + = 0.05 and n + = 0.2 are compared with the energy differences found in figure 
13 between the cases Re = 670 and 1410. All of the profiles are normalized to have the 
same slope at the wall (and therefore the same rms of the wall shear stress). The similarity 
lends support to a model based on the laminar Stokes solution with typical inactive-motion 
frequencies of the order of n + = 0.1. Bradshaw (1967) proposed a model based on small 
quasisteady perturbations to the turbulent velocity profile. This model predicts for the 
inactive motion an rms proportional to (t/^ + y+dUq /dy + ), where (y + ) is the usual 
law of the wall. This quantity is also plotted in figure 16, with the same normalization, and 
is seen to be much larger than the other curves except very near the wall. This shows that 
the hypothesis (implicit in Bradshaw’s model) that the velocity profile is in equilibrium 
with the wall shear stress is not satisfied at the present frequencies, except within a few 
wall units from the wall. The perturbed-law-of-the-wall model may supersede the laminar 
model but only at much higher Reynolds number, when sufficiently low frequencies n + 
carry energy, and then only for the u component. 

The behavior of the Reynolds stresses will now be compared with the law given by 
Townsend (1976, p. 154) and Perry et al. (1985): 

< u 2 > + = B x - Ailog(~;), < w 2 > + = B 2 - A 2 log(^), <v 2 > + =A 3 . (12) 

This law is thought to apply only in the fully turbulent region and at sufficiently high 
Reynolds numbers (as already mentioned, Perry et al. also proposed a viscous correction 
to extend this law to lower Reynolds numbers). The law implies that the stresses depend 
only on y/6, which was well verified (after viscous correction, see figure 12). At a fixed value 
of y + , a term proportional to log{du T J i>) enters the formula for < u 2 > and < w 2 > (but 
not < v 2 > or < uv >). The results in figure 13 are quite consistent with this; notice the 
much smaller rise of < v 2 >. Equation (12) predicts constant profiles in figure 16, which 
is consistent since it applies only outside the region affected by the no-slip condition. 
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Figure 17 shows a final test of (12), comparing the y-dependence of the corrected 
Re = 1410 results with the logarithmic behavior predicted by (12), using A\ = 1.1, A 2 = 
0.66, A 3 = 1.75, B\ = 2., and B 2 = 1.1. The first two values were chosen after inspection 
of the spectra (§3.4); the last three values were chosen empirically. Perry et al. (1985) 
give slightly different values: A 3 1.5 to 1.9, B\ = 2.48, B 2 = 1.12. The agreement in 
figure 17 is moderately good. Recall that the wall region, y + < 50, y/S < 0.075, should 
be excluded, and that the slope of the straight lines [A j , >42, and 0) was not adjusted to 
obtain the best agreement. Because of the insufficient Reynolds number, the present data 
cannot provide a definitive confirmation of (12), but they are consistent with it as shown 
in figures 12, 13, and 17. 

An unexpected implication of (12) is that at a given value of y+, the stresses < u 2 > + 
and < w 2 > + not only depend on the Reynolds number, but even tend to 00 like log (6 + ) 
as the Reynolds number tends to 00 (however, the log-layer overlap condition requires 
log( 6 + ) = Uoo/ilt + B , where B is a constant, so that < u 2 > /U^ and < w 2 > /U^ still 
tend to 0). Such a behavior has not been indicated by experiments, but measurements in 
the wall layer at very high Reynolds numbers are difficult and may not be accurate enough 
to reveal a slow, logarithmic divergence. Direct simulations at Reynolds numbers higher 
than the present ones have not been conducted. Presumably, as the Reynolds number 
increases, the added energy will be independent of y + away from the wall and “damped” 
by a factor similar to the ones in figure 16 near the wall. This damping will extend to higher 
values of y + as the frequency n + decreases so that the effect at a y + of 10, for instance, 
will be less than in the log layer. As a result, the shape of the Reynolds-stress profiles 
near the wall will be progressively altered; for instance the peak value of < u 2 > + will be 
displaced from its usual location into the log layer. In Klebanoff’s (1954) measurements 
at Re ss 7500, the peak value of the rms of u is about 2.98 u T and occurs at y + w 22. 
These values are significantly higher than the accepted values, which are about 2.7 and 
13. Finally, note that the diagonal elements of the Reynolds-stress tensor can tend to 00 , 
compared with the off-diagonal elements, without violating the condition of realizability. 
The correlation coefficients simply tend to 0. 

Another implication of the theory is that in (12) the constants B\ and B 2 are not 
universal, so that the Reynolds stresses at a fixed y + depend not only on <5 + , but also 
on the type of flow: e.g., boundary layer, pipe, and channel (the definition of 6 is not 
universal either). Thus (12) conflicts in many ways with the concept of a law of the wall 
for the Reynolds stresses. Bradshaw (1967) used his quasi-steady perturbation model to 
show that the mean-velocity profile has very little sensitivity to the inactive motion, which 
could explain why it does not deviate from the law of the wall. 

The size of the flow structures at different Reynolds numbers was examined in relation 
with the failure of the law of the wall for the statistical quantities. Figures 18 and 19 show 
contours of constant vorticity magnitude at Re — 300 and 1410, normalized by yJu z T /bv 
(recall that this normalization yielded a collapse of the vorticity rms). Cross-sections of the 
flow by a streamwise plane, a vertical spanwise plane and oblique planes at ±45° are taken 
in the manner of Falco (1977) and Head and Bandyopadhyay (1981). Similar figures were 
obtained from numerical simulations in channels (e.g., Deardorff 1970, Schumann 1975, 
Moin and Kim 1982.). Note the depth of the irrotational “valleys” and the sharpness of the 
irrotational-rotational interface (compared with the boundary-layer thickness), especially 
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at the higher Reynolds number. This suggests that the entrainment process is at least 
qualitatively reproduced within the multiple-scale approximation (essentially because V 
is negative: V ss -0.0251 /qo at y — 6 ). Note also the absence of numerical oscillations 
in the irrotational region, as could be caused by inadequate resolution. The figures show 
that the “typical eddies” are smaller (relative to the boundary-layer thickness) at higher 
Reynolds number, as shown by the experiments. The similarity between the present results 
at R .0 = 1410 and Head and Bandyopadhyay’s photographs at Re — 1700 is good. It is 
however difficult to ascertain that the eddy scales in figures 18 and 19 are exactly the same, 
in wall units, as Head and Bandyopadhyay’s model implies. 

3.4. Analysis of the spectra 

The Re — 1410 case is used for a more detailed study in order to minimize the low- 
Reynolds-number effects. In a wall-bounded flow, wave numbers can be nondimensional- 
ized in at least four ways; one can use the Kolmogorov length scale 1^ = (u 3 /e)^ 4 , the 
wall length scale i//u r , the distance from the wall y, or the thickness of the boundary 
layer, 6 (Perry et al. 1985). The presentation of the results will make use of the different 
normalizations and owes much to Perry et al.' s work (1985, 1986). Their theory is based on 
certain assumptions about the shape of the eddies, and on dimensional analysis. Bradshaw 
(1967) also presented a theory based on dimensional analysis. Perry et al. tentatively de- 
fine the “fully turbulent region” by y + > 100, y/6 < 0.15. Partly because of the moderate 
value of the Reynolds number, a less conservative definition was adopted here: y + > 50, 
y/6 < 0.3. The study of the spectra indicate that this definition is not unreasonable. 

Figure 20 shows spectra with Kolmogorov scaling. The maximum nondimensional wave 
numbers kl k are about 0.3 and 0.8 in the x and z directions, respectively. The values of 
y + are 100 and 200 (recall that « 660). The spectra from the twp planes collapse 
quite well for kl k larger than about 0.05. Deviations appear at the highest wave numbers, 
probably because of numerical truncation; there is even a slight turnup of the z-spectra. 
In the x direction, the spectrum of < u 2 > has a significant inertial range, about half 
a decade. For the other components, the inertial ranges are short or insignificant. The 
Kolmogorov constant was computed as the maximum of E(k)k b / 3 . As for the log layer, 
this definition has the advantage of being rigorous, but one should expect the Kolmogorov 
constants to be slightly overestimated especially if the spectra are noisy. The value for 
< u 2 > is 0.55 ± 0.05, which is in good agreement with the accepted value of about 0.5 
(Townsend 1976). The transverse Kolmogorov constants are 0.62 and 0.7, with similar 
uncertainties, for < v 2 > and < w 2 >, respectively. Recall that in isotropic turbulence, 
the transverse Kolmogorov constants are larger than the longitudinal one by a factor of 
4/3. 

The behavior of the z-spectra (figure 20b) is similar to the x-spectra for kl k larger than 
0.1, but for lower wave numbers no inertial range is found. The spectra do not remain 
tangent to the /c~ 5 / 3 law; instead they peel off for k z l k lower than about 0.05. Furthermore 
the Kolmogorov constants, defined as before, deviate from the accepted values more than 
in the x direction: they are about 1.4, 0.7, and 0.6 for < u 2 >, < v 2 >, and < w 2 >, 
respectively. It seems that a much higher Reynolds number is needed for the inertial-range 
behavior to be set in the z direction. 
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Figure 21 shows a test of the isotropy of the small scales. The one-dimensional spectra 
of an isotropic field satisfy £ 2 = £3 = ( E\ — k\dE\ / dk \) /2 where k\ is the wave number, 
E\ the spectral density in that direction, and £ 2 and £3 the densities in the transverse 
directions (Townsend 1976). These three quantities are compared in figure 21, in the z 
direction; the x-direction spectra show the same trend but are noisier. In figure 21a, with 
y + = 200, the small scales are seen to satisfy the isotropy conditions rather well for k\lk 
larger than about 0.1. In contrast, in figure 21b, with y + = 40, the spectra do not show 
isotropy. A similar test was conducted by Klebanoff (1954). The lack of isotropy of even the 
small scales near the wall is explained by the fact that the macroscales of the turbulence, 
namely y and v/u T , are not much larger than the Kolmogorov scale Ik- The size of the 
buffer-layer vortices, which are obviously nonisotropic, nonuniversal structures, is about 
15i '/u T whereas Ik is approximately equal to u/u T . This implies that the fundamental 
assumption of the large-eddy-simulation method, the existence of universal, isotropic small 
eddies, is not valid in the wall region. 

In theory, the spectra allow computation of the value of the constant C j (§3.3). As 
presented by Perry et al. (1985, 1986) the correction is C\/\/y+, (4/3)Cj/ \fy+ , and 
(4 / 3)C 1 / \/y+ for < u 2 > + , < v 2 > + , and < w 2 > + , respectively. It seems disturbing 
that this correction is not isotropic. To obtain the 4/3 factor, ones needs to restrict 
one’s attention to streamwise spectra, and to make the nontrivial assumption that the 
transverse spectral densities equal 4/3 times the longitudinal spectral density even in the 
viscous peel-off range. This assumption is not compatible with the equation £ 2 = £3 = 
(£1 - k\dEi/dki)/2 and the fact that in that range, dEi/dky < (-5/3)£i//: 1 . In fact, 
using the isotropy equation, one can easily show that the integral is the same for the three 
velocity components. For the present study (including figure 12) it was decided to discard 
the 4/3 factors. 

The values of the constants C\ and C j are not easy to establish. They are related 
by C j = C\ j yfiz since in the constant-stress layer c + = l/(/cy + ). Perry et al. quote 
C\ = 6.08, which implies C 1 = 3.89. They obtained this value by setting the spectrum to 
0 beyond a peeloff wave number klk = M, which could yield a rather crude over-estimate 
of the energy difference; M is also ill-defined. From the spectra in figure 20b, one obtains 
the estimate C x « 2. However the spectra are thought to be somewhat inaccurate (too 
high) near the numerical cutoff, and the value of the Kolmogorov constant is not firmly 
established, so that this estimate is not very reliable either. To obtain a reliable estimate 
one would probably need a simulation that exceeds the Kolmogorov wave number in all 
directions, and shows a convincing inertial range. Studies by Grant and Pao (see Hinze 
1975) suggest values of about 5 and 3.4, respectively. The value C j = 2.8±0.2 was obtained 
empirically by optimizing the collapse of the profiles in figure 12. 

In figure 22, 6 is used as a length scale. Perry et al. (1985) predict that for low 
wave numbers ( k x 6 less than about 3), the spectra from different horizontal planes will 
collapse for u and w , but not v. Although the spectra are noisy especially at low wave 
numbers, figure 22 shows encouraging agreement with the theory. The collapse is good 
for u in both directions and rather good for w. The v-spectra are higher for higher 
values of y, especially in the z direction. The agreement with the experimental spectra of 
Bradshaw (1967) and Perry et al. (1985) is good for k x 6 larger than about 1. However for 
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lower wave numbers, the computed spectra do not level off like the experimental spectra. 
This is especially apparent for the spectra of u. Each set of data has its own sources of 
error. The experimenters converted frequency spectra into spatial spectra using Taylor’s 
hypothesis, which is least valid for low wave-numbers (see Perry et al.'s discussion of the 
errors involved). The simulations directly yield spatial spectra, so these errors are avoided. 
On the other hand, the multiple-scale approximation may be showing its limitations here. 
In a true spatially-developing flow, the incoming flow has a smaller thickness and therefore 
the long waves are less energetic. Since this effect is suppressed by the approximations that 
were made one would expect the energy of the long waves to be slightly over-predicted, 
which is what the comparison between figure 22 and the experimental results suggests (the 
crude treatment of the straining effects, dicussed in §2.2, could also play a role). 

Finally figure 23 shows the spectra normalized with y as the length scale. This scaling 
provides the richest behavior in Perry et al.'s (1985) theory: a collapse of spectra from 
different values of y, with a pivot point at k x y of the order of 1, a k~ l range on the 
left (except for v which is expected to be constant), and the k x 5 ^ Z range on the right. 
Experimental results generally confirm this picture (Klebanoff 1954, Perry et al. 1985). 
In figure 23, segments with slope —1 and —5/3 are tentatively drawn on all the plots. In 
the x direction, the results show a good collapse on the right of the pivot, with the curves 
corresponding to lower values of y peeling off first. The degree of agreement with a fc~ 5 / 3 
law was discussed earlier. In the z direction the collapse, like the k~ 5 / 3 behavior, may be 
prevented by the insufficient range of scales. However, there is no reason why the theory 
should not apply to spanwise spectra. 

On the left of the pivot, the collapse and the agreement with the A; -1 model are convinc- 
ing only for < u 2 > and < — uv >, in the x direction. A conclusive check of the theory is 
made difficult by the statistical noise, the finite values of the periods A z and A z , and the 
limited range of scales due to the moderate value of the Reynolds number. From figure 23 
one can estimate the constants multiplying for the various components. For < — uv >, 
it is 0.18, which is in excellent agreement with Klebanoff’s (1954) measurements. Using 
Perry et al.'s notation, the computed values (which are probably slightly overestimated) 
are: for u, A\ = 1.1; for w, A 2 = 0.66. These constants entered (12) (in fact (12) is a 
direct consequence of the existence of a fc _x range). Perry et al. (1985, note that their 
values are switched in their Table 1) quote Ai = 1.03, A<i — 0.73. In a pipe flow Perry 
et al. (1986) measured A] — 0.90; Klebanoff (1954) obtains A\ « 0.85. These constants 
are thought to be universal. The agreement between measured and computed values is 
acceptable. 


3.5. Reynolds-stress budget equations 

The various terms in the budget of the four nonzero Reynolds stresses are shown in figure 
24, near the wall, and in figure 25 away from the wall. The contribution of the growth 
(or advection) terms is negligible near the wall (recall that V is zero in that region) and 
becomes noticeable only for y/6 larger than about 0.5. Conversely, the viscous diffusion is 
negligible for y + > 25. The terms sum up to 0 very well near the wall; in the outer region 
the residuals are less than 0.5 in the units of figure 25. As expected the near-wall behavior 
is very similar to that observed in a channel flow (Moser and Moin 1984). As the Reynolds 
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number increases, there is a weak trend for most of the quantities in figure 24 to increase. 
Furthermore there is a very noticeable enhancement (up to 20% increase from Re — 300 
to 1410) of the dissipation and of the viscous diffusion for y + smaller than about 10 (as 
part of the inactive- motion theory Bradshaw (1967) predicted that the dissipation and the 
turbulent diffusion would be enhanced). Since the inactive motion has long time scales, one 
expects it to contribute relatively less to the energy budget than to the energy itself. The 
enhancement of the viscous diffusion is of course linked to the increase of the Reynolds 
stresses, seen in figure 13. The profiles in figure 24 should be considered as representative, 
but not as universal. 

In figure 25 the tefms are normalized by uf/6. Again, one should not take the results 
as universal. At higher Reynolds numbers the dissipation tensor slowly becomes more 
isotropic. This evolution of the dissipation terms is compensated for by the pressure terms 
and the turbulent-diffusion terms, in roughly equal parts. Note that numerical truncation, 
by suppressing small-scale structures, tends to artificially increase the anisotropy of the 
dissipation tensor; the refined simulation that was done at — 500 showed slightly less 
dissipation-tensor anisotropy than the basic simulation. However the differences between 
the two sets of results were of the order of only 0.03 for the nondimensional anisotropy 
(which is defined as the ratio of the various terms to the trace of the tensor, and takes 
values of the order of 1/3). For the shear stress, the production and pressure terms 
both increase significantly with Reynolds number. Thus unless it explicitly includes low- 
Reynolds-number effects, a turbulence model should not be expected to match the present 
results very closely. For the total energy, the production and the dissipation are in balance 
over most of the boundary layer; but for separate components, the pressure term is often 
significant. For y/6 between about 0.5 and 1.2, the turbulent-diffusion and growth terms 
become important. In this flow the contribution of the growth terms to the Reynolds-stress 
budget is consistently negative. 


4. Conclusions 

Direct numerical simulations of the boundary layer on a flat plate were performed at 
four Reynolds numbers, including a value significantly larger than in any previous direct 
simulations. A systematic multiple-scale procedure was used to approximate the local 
effects of the streamwise growth of the flow. This approximation seems to be successful. 
The numerical truncation errors were estimated and shown to be much smaller than the 
relevant variations of the results. In general the agreement with experimental results was 
good. The most significant discrepancy is in the value of the computed friction coefficient 
at the highest Reynolds number, which is too large by about 5%. 

The mean-velocity and shear-stress profiles behaved as predicted by Coles’ analysis, 
except for the strength of the wake. At very low Reynolds numbers, Re less than about 
600, the logarithmic layer disappears, and the Clauser-plot method yields incorrectly high 
values of the friction velocity. A definition of the boundary-layer thickness 6 as an integral 
of the total shear stress was introduced, and produced a very good collapse of the stress 
profiles from different Reynolds numbers. The total stress appears to approach the wall 
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with a finite slope, in contrast with the parabolic behavior that is usually assumed. The 
other Reynolds stresses also collapsed very well, away from the wall, after a variant of 
Perry’s viscous correction was applied. When the bound ary- layer thickness 6 is used as a 
length scale, the energy of the small-scale motion increases with Reynolds number. 

Near the wall the Reynolds stresses < u 2 > and < w 2 >, unlike the mean velocity 
and the shear stress, did not collapse when normalized with u and u T and showed a 
strong trend to increase, consistent with Townsend and Perry’s theories. This model even 
predicts that the stresses at fixed y + tend to infinity like the logarithm of the outer-flow 
Reynolds number, and depend on the type of outer flow. When the wall length scale 
v /u T is used, the large-scale motion is responsible for the increased energy. The waves fit 
the description of “inactive motion” given by Townsend and Bradshaw. A simple model, 
based on an oscillating near-wall layer, predicted the proper trend for the stress variation 
with Reynolds number, but a convincing quantitative extension of the model into the 
wall region remains to be found. The spectra showed encouraging agreement with various 
theories, including Kolmogorov’s and Perry’s, and yielded satisfactory values for some of 
the universal constants. Extensive data were provided for the development of turbulence 
models both near the wall and in the outer region, and the Reynolds- number dependence 
of the major terms was discussed. The anisotropy of the dissipation tensor was found to 
be significant even at the highest Reynolds number reached. 


The author had useful discussions with Drs. J. Kim, N. Mansour, P. Moin, R. Rogallo 
and A. Wray (NASA Ames Research Center), and with Prof. P. Bradshaw (Imperial 
College, London). Dr. Mansour reviewed the manuscript. 
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APPENDIX 


We examine the classical two-layer model of the flat-plate turbulent boundary layer 
(Townsend 1956, p. 229), especially the local momentum equation. The assumptions are, 
near the wall: 

U + = f(y + ), (Ala) 

t + = 1, (Alb) 

and away from the wall: 

(U - Uoo) + = g(y/ 6 ), (A2a) 

t + = h(y/ 6 ), (A2b) 

for some nondimensional functions /, g and h. The boundary-layer momentum and con- 
tinuity equations are 


UU x + VUy= — , U x + V y = 0. (A3) 

dy 

In the wall region, which includes the log layer, the momentum equation becomes (after 
using some algebra and continuity): 


h' = 6 -^f\ (A4) 

Strictly speaking this contradicts (A2b), because the right-hand-side of (A4) is not just 
a function of y/S. Is it compatible with a parabolic behavior of g, in particular as the 
Reynolds number tends to oo? When this happens 6 u TX /u T tends to 0 and / 2 tends to oo 
(at fixed y/ 6 ); we therefore need better estimates. We shall use the overlap condition and 
the total-momentum equation. 

Equations (Ala) and (A2a) overlap in the log layer. By differentiating the overlap 
condition with respect to x, we obtain (Coles 1956) 

^ + ^(l + «^) = 0. (A5) 

6 u T u T 

The momentum equation is dO /dx = u^/U^. Using (Ala) and (A2a) we obtain 


e = d x £-6 + d 2 ^s + + D. 


Ur. 


u 2 

oo 


■f/n 


uu T 

[ Tp 

oo 


(A6) 


for some constants D\ « 3.9, D 2 ~ -24, D 3 « 52, and D 4 . Clauser (1954) gives the value 
3.6 for D\ (A in his notation). 

We compute dO j dx, neglecting the last term in (A6) which is very small, and using (A5) 
to eliminate 6 X : 

SUtx 
U r 


-kD\ + (!-«:) 


D2U t 


K 

U 2 

^ OO 


(A7) 
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so that (A4) becomes 


1 


ti = 


(*£>,- (!-*)%*■) 



(A8) 


We consider a fixed value of y/6 (but within the log layer) and let the Reynolds number 
tend to oo. The ratio u T /U oq tends to 0 and l/’/f/ 0 o tends to 1. Therefore the leading term 
is —1/(kDi) (« -0.6) and is a constant. This means that up to the edge of the log layer 
r + (y) has a flat part with finite slope —0.6/6. Townsend (p. 255, 1956 edition) obtains 
the same behavior, with slope about -0.5. Curiously, he removed this result for the 1976 
edition. At moderate Reynolds numbers, both factors in (A8) get smaller. Klebanoff’s 
1954 results and the present results show about -0.5 for h' up to y/6 « 0.2. 
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Figure 1. Sketch of a spatially-developing boundary layer, a) geometry, 
vectors; — streamlines; coordinate lines, b) turbulent signal. 
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Figure 2. Resolution check at Rg* = 500. coarse resolution; — all other cases. 
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Figure 3. Streamwise thickening of the boundary layer, o displacement thickness 6~; 
• momentum thickness 0; slope c//2. 
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Figure 6. Strength of the wake, A U + . Present results: • based on apparent log law; 
o based on normal log law. — Coles; + other experiments. 



Figure 7. Shear-stress distribution. — R$ = 300; Re — 670; Re — 1410. 

Upper curve, total stress; lower curve, Reynolds stress. • Re « 7500 (Klebanoff 1954). 
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Figure 9. Stresses and production near the wall. — Re — 300; R$ = 670; 

Re — 1410. ... l//cy + . A Reynolds stress; x total stress; o production. • Kim et al. 

1968. 



Figure 10. Turbulence intensities. Upper curve u rma /u T ; lower curve v rm3 /u T ; middle 
curve w rms /u T . — Re — 300; R$ = 670; Re = 1410. • Re « 7500 (Klebanoff 







Figure 12. Turbulence intensities, corrected to infinite Reynolds number. Upper curve 

Urms / ^r! lower curve t/ rm j/«ri middle curve w rms /u T . — R$ — 300; R$ = 670; 

Re = 1410. • R$ « 7500 (Klebanoff 1954). 



Figure 13. Turbulence intensities. Upper curve u rma /u r ; lower curve v rms /u T ; middle 
curve tv rrns /u T . — Re = 300; Re = 670; Re — 1410. 
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Figure 16. Profiles of inactive-motion rms. Energy difference between R e = 670 and 

10: — < u 2 >, < w 2 >. Stokes-layer model: + n + = 0.2, x n+ = 0.05. 

- law-of-the-wall model. 


Figure 17. Reynolds stresses, extrapolated from Re = 1410 to oo 
— - w +2 . • computed values; A equation (11). 
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Figure 24. Concluded. 
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Figure 25. Reynolds-stress budget away from the wall, a) u z \ b) v , c) w ; 

d) -uv; e) u 2 + v 2 + w 2 . -o- production; — turbulent diffusion; growth terms 

- - — dissipation; ... pressure. 
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